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Abstract 

The critical behaviour of the randomly spin-diluted Ising model in 
two space dimensions is investigated by a new method which com- 
bines a grand ensemble approach to disordered systems proposed by 
Morita with the phenomenological renormalization group scheme of 
Nightingale. Accurate approximations for the phase diagram and for 
the connectivity length exponent of the percolation transition are ob- 
tained. Our results suggest that the thermal phase transition of the 
disordered system might be different from that of the pure system: we 
observe a continuous variation of critical exponents with the density p 
of magnetic impurities, respecting, however, weak universality in the 
sense that r\ and 7/1/ do not depend on p while 7 and v separately 
do. Our results are in qualitative and quantitative agreement with a 
recent Monte-Carlo study. 
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The present contribution is concerned with the critical behaviour of the 
randomly spin-diluted (SD) Ising model in two space dimensions. For this 
model, the Harris criterion [1] according to which small amounts of disorder 
do (do not) change the nature of the phase transition, if the corresponding 
pure system's specific heat exponent is positive (negative), is inconclusive, 
since a pure = 0. The critical behaviour of the disordered system has been 
a subject of debate for many years; see e.g. [2] and references therein, and 
[3-10]. 

Currently, there appears to be widespread consensus — based upon a 
field-theoretic reformulation of the system with weak bond-disorder (BD) in 
the critical regime — that, in this limit, critical exponents of the disordered 
system are the same as those of the pure system, albeit modified by loga- 
rithmic corrections [3-6]. Numerical simulations of BD systems [7,8] indicate 
that such modifications through logarithmic corrections would persist even 
in the strongly disordered regime. This consensus is, on the other hand, 
being questioned in [9], where a more complicated phase transition, and a 
non-divergent specific heat are predicted. Lastly, a very recent large-scale 
Monte-Carlo study of the SD Ising system [10] produces results which are at 
variance with the findings of [3-9], in that it gives critical exponents which 
clearly vary with the density p of occupied sites. 

It is this investigation, in particular, which finally convinced us to put 
our own transfer-matrix analysis of the SD system — both the method, and 
its results — to public discussion, even though, or rather because we are 
aware that a deeper understanding of our approach [11] would certainly still 
be welcome. In essence, we find that there is full qualitative and quantitative 
agreement between the Monte-Carlo results of Kim and Patrascioiu [10] and 
ours. 

We study the SD system by a new method [11] which combines Morita's 
grand ensemble approach to disordered systems [12] with phenomenological 
renormalization [13]. We begin by briefly describing Morita's method of 
configuration averaging in systems with quenched randomness. Within this 
approach, a hierarchy of approximating equilibrium systems is constructed, 
whose critical behaviour in turn is obtained by phenomenological renormal- 
ization. We then state the main results, relegating a detailed description of 
the method to a separate publication. 

Consider a disordered system described by the Hamiltonian H(ct\k), 
where k denotes the quenched disorder configuration, and a the set of 
dynamic variables of the system. Morita's approach avoids configuration 
averaging of the free energy by working in an enlarged phase space in 
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which k, too, is a dynamic variable, and by introducing a potential (ft(n) 
chosen such that the new system with hamiltonian H^(<t,k) = H(a\n) + 
<P(k) exhibits thermodynamic equilibrium properties identical to the non- 
equilibrium properties of the original quenched system. To achieve this, the 
distribution p^(a, k) generated by H^(a, k, ) must be constructed such that 
it satisfies 

p*(<7, «) = ^ ex P [-/?^(a, «)] = ||| exp[-/3i7(a|K)] (1) 

for all (a, k), or equivalently —[34>{k.) = ln[g(/t)/Z(K)] + \nZ^. Here and 
are partition functions of the grand ensemble and the quenched system 
with fixed impurity configuration k, respectively; g(/c) is the probability 
distribution describing the quenched disorder. If (ft is normalized such that 
its configuration average vanishes [14], Eq. (1) implies 

lnZ^ =< In Z(k) > q - < \nq(n) > q , (2) 

where < .. > q denotes an average over the quenched disorder. That is, lnZ^ 
gives the Brout free energy plus an irrelevant contribution of an entropy of 
mixing. 

Eqs. (1) and (2) are not very useful for practical calculations. In order 
to utilize Morita's ideas, one has to find a representation for 4>(k) which 
is adapted to the problem at hand, and one will usually have to resort 
to approximations [11,15-17]. For the spin-diluted Ising model, one may 
expand according to 

P<t>(n) = A + Ai ^ fci + A 2 ^ fcifcj + . . . + A P ^ IJ fci + . . . (3) 

i p ieP 

where the ki are occupation numbers taking the value 1 or if in k the site 
i is occupied or empty [18]. The first sum in (3) is over all lattice sites, 
the second over nearest neighbour pairs, and the third over all elementary 
plaquettes of the system. Each term serves to control one moment of the 
probability distribution p^(k) = J2<jP^( a i K ) " i- e - an expectation under p^ 
of some product of occupation numbers. The couplings Ai, A2, ... have to be 
determined as functions of temperature and field such that the moments of 
P^(k) coincide with those of q(n). The constant Ao is required to achieve 
< <t>{K) > q = 0. 

For a randomly spin-diluted system with density p of magnetic impuri- 
ties, one would have to determine the couplings of (ft so that 

< ki >0= p , < kikj >0= p 2 <]Jki >0= p |p| , (4) 
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and so on [11,16]. Here < .. >0 denotes an average with respect to p^(k) 
and |P| the size of the elementary plaquette, which for the square lattice is 
4. Eqs. (4) constitute an infinite set of equations for the couplings Ai, A2, ... 
of the potential <j>. To obtain a full solution is, in general, impossible. In- 
terpreting (4) as a set of constraints imposed on the thermal motion of the 
magnetic impurities, one may, however, set up a systematic scheme of ap- 
proximations by letting only finite subsets of this set of constraints become 
operative [11,16]. Implementing only < ki >0= p , one would describe an 
annealed system at density p of magnetic impurities. Such a system would 
condense and order at low temperatures and would thus provide a rather 
poor description of the quenched system. If, in addition, one fixes nearest 
neighbour correlations < kikj >^ of particle locations at their quenched 
value p 2 , the condensation phenomenon no longer occurs and the system 
exhibits, e.g., a percolation transition. Such a system is already a serious 
candidate for the description of fully frozem-in disorder. In this manner, one 
arrives at increasingly accurate descriptions of the quenched system as more 
and more constraints are taken into account, until eventually one would ob- 
tain an exact description of the original disordered system. The hope, of 
course, is that already rather simple approximations in this hierarchy might 
belong to the "univerality class" of the quenched system. 

Except for the one dimensional system [11], this program cannot be car- 
ried through to the end, and for d > 2 we know of no exact solution of even 
the simplest approximating systems. At all levels of approximation, though, 
one has to do with translationally invariant equilibrium systems. Their crit- 
ical behaviour could, for instance, be obtained by standard RG methods, 
were it not for the fact that the couplings of <p are only determined through 
a set of constraints. Since it is far from obvious how these constraints 
should be transformed under rescaling (see, however, [16]), we decided to 
use phenomenological renormalization [13], which avoids this problem alto- 
gether. No explicit RG transformation in the space of couplings need be 
constructed. Given p, one just solves the system in strip geometries - with 
the appropriate set of constraints (4) imposed [11]. 

If £m(p, T) denotes the correlation length of an (M x oo)-system at den- 
sity p, temperature T, and field H = 0, an approximation to the critical 
point T c is given by the fixed point of the phenomenological RG relation 
£m(p,T)/M = £m'(p, T')/M'. An approximation to the correlation length 
exponent v may be obtained by linearizing this relation about its fixed point 
[13]. Note that p is a parameter of this procedure. By varying p, we can 
obtain the phase boundary T c (p), and the thermal correlation length expo- 
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nent v t along the critical line. Analogous relations can be used to study the 
percolation transition at T = 0, p = p c . 

We have studied four approximating systems, named (a) - (d). In system 
(a), only the first and second constraint displayed in (4) are imposed on the 
system. Due to the anisotropy of the strip geometry, correlations parallel and 
perpendicular to the strip turn out to be different. We have thus introduced 
system (b) where correlations parallel and perpendicular to the strip are 
treated as separate constraints. Systems (c) and (d) are obtained from (a) 
and (b) by fixing, in addition, correlations around each elementary plaquette 
at p 4 . Note that the results obtained from systems (a) and (b) or (c) and 
(d) should approach each other as the strip widths M and M' go to infinity. 
This can serve as a valuable consistency check for extrapolations of critical 
parameters obtained from fixed points of the phenomenological RG relations 
to the infinite system values. 

We now turn to the results. We have computed critical temperatures 
T c {p) for various densities p < 1. The values obtained for systems (a) - (d) 
agree very well down to p ~ 0.75. Below p ~ 0.75, critical temperatures of 
the approximating systems with the plaquette-constraint imposed turn out 
to be slightly larger than those of systems (a) and (b) without this constraint. 
Extrapolating a parabolic interpolation of the T c {p) values for p = 0.9, 
p = 0.95, and p = 0.975 to p = 1, we obtain T c (p) -> 2.2674, as p -> 1, 
which is off the mark by less than 0.1%, and T c (p) _1 dT c (p) / dp — > 1.579, as 
p — > 1, which differs from the exact result 1.535 [19] by less than 3%. We 
have also checked the scaling form of the transition line near the percolation 
transition at p c , exp[— 2 J/kBT c (p)] ~ (p — p c )^, with (p = v t (p c )/v p , and 
we find if = 1 to within less than 2% for all four approximating systems 
[11]. Here u p denotes the connectivity-length exponent of the percolation 
transition. The percolation threshold p c itself (p c ~ 0.593 [20]) is correctly 
reproduced to within 3% and 1% by systems (a),(b) and (c),(d) respectively; 
see table I. 

Fig. 1 shows our results for the correlation length exponent v of system 
(a), based on extrapolations from strip widths up to M = 8 and M = 10 
for the thermal and the percolation transitions, respectively. They clearly 
show a variation with the density p of magnetic impurities. For p = 0.7 and 
p = 0.9, table I shows that this variation persists, as we impose further con- 
straints in approximating systems (b)-(d). Note that the correlation length 
exponents are always mutually consistent for these four approximating sys- 
tems. For the percolation transition, they agree well with the expected exact 
result i/p = 4/3 of den Nijs [26]. Also, the isotropic and anisotropic vari- 
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ants in the two groups of systems, (a),(b) and (c),(d), always give consistent 
results (as they should). While all four systems agree in "universal" char- 
acteristics of the phase transition, u p or u t , the first group can differ from 
the second as far as non-universal parameters of the phase transitions, i.e. 
p c or T c (p), are concerned. 

Table 1: Extrapolated critical parameters for the percolation transition, and the 
thermal phase transition at densities p = 0.7 and p = 0.9 of magnetic impurities. 
Numbers in brackets give the estimated error of the last displayed digit of the 
preceding quantity. 





percolation 


P = 


0.7 


P = 


0.9 




Pc 


u p 


T 

± c 


v t 


T 

± c 


ft 


(a) 


0.609(1) 


1.33(2) 


1.051(1) 


1.31(2) 


1.902(1) 


1.13(1) 


(b) 


0.609(1) 


1.33(2) 


1.050(1) 


1.30(2) 


1.901(1) 


1.13(1) 


(c) 


0.587(1) 


1.34(2) 


1.081(1) 


1.29(2) 


1.901(1) 


1.12(1) 


(d) 


0.587(1) 


1.33(2) 


1.080(1) 


1.30(2) 


1.901(1) 


1.13(1) 



The variation of critical exponents with the density p is at variance with 
what would be expected from universality arguments, and one might well 
argue it to be an artefact of our simple approximations. So independent 
evidence is clearly welcome. Therefore, we turned to standard real-space 
renormalization group (RSRG) calculations for this system [2], and noted 
that they, too, are able to produce similar results: Some time ago, Tsal- 
lis and Levy [21] used the RSRG to compute phase diagrams of bond di- 
luted models, and they propose, inter alia, a set of RSRG transformations 
parametrized by the density p of bonds. Applying this idea to the SD prob- 
lem, and using techniques of [22], we obtain the full curves shown in Fig. 1. 
While such an approach may still seem ad hoc, we find that, last but not 
least, very recent large scale Monte-Carlo simulations of the SD system [10] 
also give a non-universal u t (p) in quantitative accord with our results, as 
can be seen in Fig. 1. 

The agreement between the results of Kim and Patrascioiu [10] and ours 
goes, in fact, much deeper: We computed susceptibilities and used the finite- 
size scaling relation xm(T c ) ~ M^/ u to determine j/u. We find j/u ~ 1.75 
as in the pure system, and independently of p, which implies that 7 itself is 
again non-universal. The amplitude Aq of the critical finite-size correlation 
length, £m(T c ) — AqM, is known to be related to the critical exponent n 
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according to Aq = l/nr] [23]. We have used this relation to determine tj, 
and we find 77 ~ 0.25 independently of p for the thermal transition, which is 
of course consistent with the p-independence of 7/1/ via the Fisher relation 
7/1/ = 2 — rj. We find this relation satisfied by all our extrapolated values for 
all our approximating systems, usually to within 1% or better. The same 
kind of weak universality [24] is also observed in [10], and it is of the same 
type as that encountered in Baxter's 8-vertex model [25]. 

At this point let us mention that some time ago Derrida et al. [27] found 
analogous non-universal behaviour of Binder's cumulant ratio in a family of 
self-dual bond disordered Ising models, which they ascribed to logarithmic 
corrections to finite—size scaling. However, as pointed out by Cardy [28], 
in order to disentangle logarithmic from power-law corrections to finite-size 
scaling, one might have to go to rather large strip widths, and our data, as 
yet, do not support the conclusion that the observed non-universality is only 
apparent and caused by logarithmic corrections to finite-size scaling; see 
Fig. 2 which shows the variation of the phenomenological critical exponent 
vm,m-i of system (a) with system size M, and the insert which adresses 
this particular point. 

The complete set of constraints (4) imposes the condition < Y\ ieu) h >^= 
pl w l for all subsets ui of the lattice, implying that occupancy of lattice sites 
is uncorrelated, or in other words that the correlation length describing the 
fcj-correlation functions vanishes. Our simple approximations only fix some 
of these correlations. Others can of course be computed and are found to 
vary only slightly with temperature and field. We have checked that the cor- 
responding correlation lengths remain small (at most a few lattice spacings) 
so that the approximating systems may be regarded as virtually quenched. 
Moreover, it can be shown [11] that in zero external field an infinite set 
of couplings of the exact potential <p vanishes, and is thus correctly taken 
into account already at the level of our simple approximating systems (It 
implies, for instance, that in zero field our system (a) already provides an 
exact description of the quenched system in 1-d). 

In summary, we have studied the critical behaviour of the spin-diluted 
Ising model by a new method which combines a grand ensemble approach 
to disordered systems with phenomenological renormalization. Where a 
comparison with exact results (concerning the phase diagram but also uni- 
versal aspects of the percolation transition) is possible, our data compare 
favourably with earlier RG analyses of the system (see e.g. [2] and Refs. 
therein). We know of no a-priori reason why the present method should be 
intrinsically better in describing the percolation transition, where we get a 
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value of v p consistent with the exact result of den Nijs, than in describing 
the thermal phase transition at p c < p < 1, where our results indicate a 
continuous variation of critical exponents with p — in a manner that re- 
spects weak universality. Indeed, a recent Monte-Carlo study [10] produces 
results completely in accord with ours, both qualitatively and quantitatively 
Hence, the the observed variation of critical exponents is most likely not an 
artefact of our grand ensemble description of quenched disorder. 

A deeper understanding of our method would certainly still be welcome. 
Let us mention in closing that it is easily adapted to the study of correlated 
disorder. Preliminary results for a 1-d system, where comparison with an 
exact solution is possible [29] do look encouraging. We are currently also 
using our method to study the wetting transition in the presence of sur- 
face disorder, which has recently been a subject of some controversy in the 
literature. 
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Figure captions 

Fig. 1: Correlation length exponent v of system (a) for various densities p 
(open squares). The leftmost open square was determined as a connectivity 
length exponent of the percolation transition. The exact result v p = 4/3 is 
displayed as a diamond at p c ~ 0.593. Results from a RSRG calculation, 
using the b=2 and b=3 decimation transformations of Yeomans and Stinch- 
combe are given as full curves. Full circles show the Monte-Carlo data of 
[10]; except near the percolation transition, where a precise location of the 
phase boundary is difficult, the agreement with our model (a) results is ex- 
cellent. 

Fig. 2: Variation of the phenomeneological critical exponent vm,m-i with 
strip width M, for approximating system (a); open squares: percolation; 
diamonds: p = 0.7, full squares: p = 0.9, open circles: pure system. The 
insert checks for the possibility of logarithmic corrections to the pure systems 
critical behaviour: The quantity 1/Am = [vm,m-i(p) — ^m,m-i(1)] 1 is 
plotted vs. ln(M). For p = 0.9, 1/Am is not monotonically increasing 
with M and levels off for large M, which is evidence against logarithmic 
corrections. For p = 0.7 the data themselves are not as conclusive as for 
p = 0.9. However, 1/Am appears to increase slower than linearly with 
ln(M) and there is a trend for the curve to level off. 
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